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ABSTRACT 

We perform detailed comparison of the semi-analytic halo model predictions 
with measurements in numerical simulations of the two and three point correla- 
tion functions (3PCF), as well as power spectrum and bispectrum. We discuss 
the accuracy and self-consistency of the halo model description of gravitational 
clustering in the non-linear regime and constrain halo model parameters. We 
exploit the recently proposed multipole expansion of three point statistics that 
expresses rotation invariance in the most natural way. This not only offers techni- 
cal advantages by reducing the integrals required for the halo model predictions, 
but amounts to a convenient way of compressing the information contained in 
the 3PCF. We find that, with an appropriate choice of the halo boundary and 
mass function cut-off, halo model predictions are in good agreement with the 
bispectrum measured in numerical simulations. However, the halo model pre- 
dicts less than the observed configuration dependence of the 3PCF on ~ Mpc 
scales. This effect is mainly due to quadrupole moment deficit, possibly related 
to the assumption of spherical halo geometry. Our analysis shows that using 
its harmonic decomposition, the full configuration dependence of the 3PCF in 
the non-linear regime can be compressed into just a few numbers, the lowest 
multipoles. Moreover, these multipoles are closely related to the highest signal 
to noise eigenmodes of the 3PCF. Therefore this estimator may simplify future 
analyses aimed at constraining cosmological and halo model parameters from 
observational data. 

Subject headings: cosmology: theory — large-scale structure of universe — meth- 
ods: statistical — numerical 



institute for Astronomy, University of Hawaii, 2680 Woodlawn Dr, Honolulu, HI 96822, USA 

2 Instituto de Ciencias del Espacio (IEEC/CSIC), Facultat de Ciencies UAB, Torre C5 -par- 2a planta, 
08193 Bellaterra (Cerdanyola), SPAIN 

3 School of Physics and Astronomy, University of Nottingham, Nottingham NG7 2RD, UK 



- 2 - 



1. Introduction 

Galaxy formation and evolution still lacks a compelling explanation from first princi- 
ples. In the absence of a successful theory, galaxy clustering can be described assuming that 
galaxies are biased tracers of the underlying dark matter distribution. Recently, the halo 
model for gravitational clustering (Neyman & Scott 1952; Peebles 1974; McClelland & Silk 
1977a,b, 1978) has been revived in cosmology as an attempt to provide an accurate picture 
of gravitational clustering in the non-linear regime as seen by high-resolution N-body simu- 
lations (Scherrer & Bertschinger 1991; Sheth & Jain 1997; Seljak 2000; Ma & Fry 2000a,b; 
Scoccimarro et al. 2001; Cooray & Hu 2001). 

Halo model provides a simple prescription for the analytic computation of N-point cor- 
relation functions which are the most widely used statistics for gravitational clustering (see 
Cooray & Sheth (2002) and references therein). In particular, Seljak (2000) suggested 
that the power spectrum of dark matter and galaxies are consistent with this approach, 
whereas Peacock & Smith (2000) reached similar conclusions using mock catalogs. Ma & 
Fry (2000a,b) derived predictions for the power spectrum, bispectrum and its Fourier coun- 
terparts, the two and three-point correlation functions, that were found to be in agreement 
with numerical simulations. Scoccimarro et al. (2001) (hereafter SSHJ), discussed the halo 
model predictions for the power spectrum, and bispectrum, as well as the higher-order cumu- 
lants, S p , and concluded that dark matter clustering in the halo model is in good agreement 
with measurements from N-body simulations on large scales (r > IMpc/h). However, the 
agreement was at the 20 % level only on smaller scales (r ^ IMpc/h). They also showed a 
first comparison with the variance from APM galaxies. 

Progress on numerical simulations has triggered a number of theoretical developments in 
the halo model. New elements such as substructure (see e.g., Sheth 2003; Sheth & Jain 2003; 
Dolney et al. 2004) and halo geometry (Jing & Suto 2002) have been recently incorporated 
into the model to provide a more realistic description of dark matter clustering. Other 
developments have focused on describing galaxy bias using the so-called halo occupation 
distribution, a prescription of filling halos with galaxies in a stochastic fashion (Berlind & 
Weinberg 2002; Berlind et al. 2003; Zehavi et al. 2004). 

Recently, Takada & Jain (2003a) have thoroughly explored halo model predictions for 
the three-point correlation function (3PCF hereafter) of dark matter and galaxies. They 
concluded that halo model predictions for the dark matter described adequately only cer- 
tain triangle configurations. They fail to match simulations for triangles with side length 
r Ki IMpc/h. Wang et al. (2004) investigated the two- and three-point functions and com- 
pared analytic results to dark matter simulations and galaxy clustering observations from 
the 2dF. They claim their predictions closely reproduce numerical results and observations. 
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In this paper we focus on the dark matter three-point correlation functions in real and 
Fourier space. In particular, our principal aim is to carry out a detailed comparison of halo 
model predictions with measurements in high resolution simulations to test the validity of 
the model and constrain halo model parameters. 

The present analysis differs from previous work in the literature in several ways: 

• We present new technology to perform the necessary integrals using multipole ex- 
pansion, in particular the dimension of the integrals is significantly reduced in our 
prescription. 

• For the first time, we test the validity of the halo model against N-body simulations 
simultaneously in real and Fourier space, at the two and three-point level. 

• We constrain halo parameters in a self-consistent fashion: all statistics are predicted 
from halo model parameters. In particular, we do not use the Smith et al. (2003) as 
an ingredient of the model, which has become common usage when constraining bias 
from observational data (see e.g., Zehavi et al. 2004; Wang et al. 2004). We show that 
using the halo model in a self-consistent way provides a good fit to N-body simulations, 
provided one leaves the halo boundary as a free parameter, and large-mass haloes, not 
present in simulations, are removed from the mass function accordingly. 

• We extend our analysis for the recently proposed multipole expansion of three-point 
statistics, which expresses rotation invariance naturally. This allows the compression 
three-point statistics into a few multipoles even in the non-linear regime. This is 
especially convenient for constraining non-Gaussianity from gravitational clustering. 

2. Halo Model 

According to the halo model picture, the non-linearly evolved dark matter distribution is 
described in terms of the clustering properties dark haloes, quasi-equilibrium objects formed 
by gravitational collapse, (see Cooray & Sheth 2002, for a review and references therein). In 
this context, correlation functions are decomposed into contributions arising from correla- 
tions among particles inhabiting dark-matter haloes. For convenience throughout this paper 
we shall closely follow SSHJ to describe the statistics of gravitational clustering within the 
halo model. 
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2.1. Halo profile 

Numerical simulations suggest that dark haloes have a universal profile (Navarro et al. 
1997; Moore et al. 1998). Here we shall adopt the NFW profile, 

P ^ = 4tt (^) cr/R v {l + cr/i2„) ^ 

and its Fourier transform, p{k) = 4irr 2 dr p(r)j (kr) , where jo(x) = sinx/x is the zeroth- 
order spherical Bessel function, takes the form (see SSHJ), 



p(k,y) = f 



sin 77 {Si //(l + c)] - Si(//)j + cos?? {C%(1 + c)] - Ci(/,)} SUl(//( ; 



r/(l + c) 



(2) 



with / = 1/ [ln(l + c) — c/ (1 + c)], T] — kR v /c, c is the concentration parameter, and Si and 
Ci are the sine and cosine integral functions. The virial radius, R v , is the characteristic 
scale which separates the inner r -1 from the outer r -3 behavior of the profile. Above we 
have assumed that the halo profile is truncated at R v , an assumption that shall be relaxed 
in order to match clustering measurements in N-body simulations (see §4). The concen- 
tration parameter is not well constrained from simulations and we adopt the standard halo 
parametrization with halo mass 

f M\~ p 

c(M) = co — (3) 



sets the non-linearity scale cr(M#) = 5 C , being a the linear rms mass fluctuation and 
S c = 1.686 is the linear overdensity required for spherical collapse. We shall assume the 
values suggested by numerical simulations cq — 9, j3 — 0.13 Bullock et al. (2001), but have 
checked that changing halo concentration parameters by as much as 50 % does not affect 
significantly our predictions within the range of scales probed by our simulations. 



2.2. Halo mass function 

A basic ingredient of the halo model is the mass function, n(M), which describes how 
many objects of mass in the range M and M+dM end up collapsing to form bound structures. 
Here we use the extension of the Press-Schechter formalism provided by Sheth & Tormen 
(1999) that accurately describes N-body results (Jenkins et al. 2001): 

n(M)MdM = p^n(u) = p^-A 1X f^- (l + g( v )-*) exp(-<y(i/)/2) (4) 
y y V 2tt 
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with 7 = din a 2 / din R, g(v) = au 2 , a = 0.707, v = 5 c /a, A = 0.322, p = 0.3, y = (R/R*) = 
(M/M,) 1 / 3 , R is the la grangian radius, and R* = R v A is the non-linear (lagrangian) scale, 
where A = 200,340 for VL m = 1,0.3 cosmologies, respectively. 



2.3. Halo clustering 

Halo correlations can be included in the model by assuming a biasing prescription be- 
tween the halo and the underlying mass distribution. We follow the approach introduced 
by Mo & White (1996) in the context of the spherical collapse model, and adopt the fitting 
formula given by Sheth & Tormen (1999): 

Kv) = 1 + sVt±+ ?P (5) 
5 C 5 c (l + g(u)p) 

and we neglect quadratic and higher-order biasing terms. We note that a non-vanishing 
quadratic bias only affects the three-point statistics (or higher orders) on large scales and 
we have checked that its contribution is a few per cent at most. 



3. Halo Model Statistics 

Halo model statistics are particularly simple to formulate in Fourier space. Complica- 
tions arising from convolutions of halo profiles become simple products in transform space. 

According to the halo model approach, the non-linear power spectrum of the mass 
fluctuations is the result of two independent contributions: one coming from the single 
halo profile auto-correlation that dominates on small-scales and another one given by the 
correlation among dark matter particles in different halos that accounts for the large-scale 
clustering. This way we can write, 

P(k) = \I ll (k)} 2 P L (k) + I 02 (k,k) (6) 

where, 

Iijih, •••,%) = J ^-n(y)bi(y) \p(k u y) ■ ■ y)\ (jr^j ( 7 ) 

where b = 1, 61 = b(u), and 6j = for i > 1 as we neglect quadratic and higher-order biasing 
terms, and Ph{k) is the linear power spectrum, for which we assume the Bond & Efstathiou 
(1984) parametrization. The two-point correlation function follows by Fourier transforming 
P(k), 

£(r) = J ^dk P(k) j (kr) (8) 
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where jo is the zero-order spherical Bessel function. 

The bispectrum, P 123 = B(k 1: fc 2 , fc 3 ), can be expressed as a sum of three-point correla- 
tions among mass particles residing in one, two, and three haloes, respectively: 



-B123 — 



AS + {hi{ki)h2{k 2 , h)P L (h) + perm(l, 2, 3)} + / 03 (fci, fc 2 , fc 3 ) (9) 



where P PT denotes the bispectrum from second order Perturbation Theory (PT). By design, 
halo model recovers weakly non-linear theory predictions on large scales, since fc->0,/n«l 
and /12, 102, 103 ~ 0, or equivalently, P P l , P i23 ~ B PT . 



The 3PCF is a triple Fourier transform of the bispectrum, 

r 3 

C(ri,r 2 ,r 3 ) = / Y[d% 5 i23 e^ kl - ri+k2 - r2+k3 - r3 )5 Dirac (k 1 + k 2 + k 3 
J i=i 



(10) 



In the context of halo models, C is a complex object: Takada & Jain (2003a) have shown 
that the evaluation of the two- and three-halo terms in real space involves eight- and twelve- 
dimensional integrals, respectively; this is untractable. They reduce the dimensionality of 
the integrals to two by switching to Fourier-space and taking a number of approximations 
in the corresponding kernels. 

Our approach is significantly different in that, although we also use Fourier-space for- 
malism, we further decompose correlation functions into harmonic multipoles. This greatly 
simplifies the exact analytic expressions, therefore we do not need to take the approximations 
used by Takada & Jain (2003a) for the 2- and 3-halo terms. We note that, for the 1-halo 
term, their estimator is expressed as a four dimensional integration and it is exact. The only 
approximation we take is that we use a finite number of multipoles. We will see that this is 
an excellent approximation, except for degenerate (isosceles) triangles, where the amplitude 
of the harmonic coefficients decrease slowly with multipole order. In our formalism, the 
3PCF can be expressed as a two-dimensional Hankel transform of the bispectrum, for each 
harmonic multipole (Szapudi 2004): 

C(ri,r 2 ,0) =X;^j^6(ri,r 2 ) P,(cos0) (11) 

where 9 is the angle between r x and r 2 , Pi is the legendre polynomial of order £, (e(r 1 ,r 2 ) is 
the harmonic 3PCF, 

Un,r 2 ) = J dh dk 2 ^^(-l) e B e (k u k 2 )j e (k iri )j e (k 2 r 2 ) (12) 
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and Be is the Legendre transform of the bispectrum, 

B e (k 1 ,k 2 ) = 2n j d(cos9) B(h,k 2 ,9) P e (cos9) . (13) 

High multipoles {£ ^> 1), only important for the isosceles triangles, can be easily computed 
employing the large-£ limit of spherical Bessel functions, je(x) m a/ tv] (2£ + 1) bmraJS- + 
1/2 -x), 

c^,r 2 ) = fAy 3 _i_ ( _ 1) ^^ ) M (14) 



\2n J (rxr 2 ) 3 \ r i r 2 

where L = £+1/2. It is convenient to define the reduced (or hierarchical) 3PCF, 

n( \ _ C{ri,r 2 ,r 3 ) 

^ e(ri)e(r 2 ) + e(r 2 )e(r 3 ) + tfr^n) 1 J 

and similarly, the reduced bispectrum in Fourier space, 

-^123 

g(A;i ' h ' h) = p(h)P(k 2 ) + p(k 2 )p(k 3 ) + p(k 3 )p( kl y (16) 

PT and the stable clustering hypothesis predict them to be weakly dependent of scale in 
the quasi-linear and non-linear regimes, respectively. Similar to eq(13), it is convenient to 
decompose the reduced 3PCF, Q(ri, r 2 , r 3 ), in its harmonic multipoles by taking its Legendre 
transform, 

Q t {r u r 2 ) = 2n J d(cos9) Q{n,r 2 ,9) P e (cos9) (17) 

where 9 is the angle between ri and r 2 . Analogous expressions can be written for the 
harmonic multipoles of the reduced bispectrum. 

Since the bispectrum itself is obtained by performing a one-dimensional integration of 
the halo correlations over the mass function, in order to compute the three-point function 
multipoles Q, as given by eq(12), four-dimensional integrals are required for each £ 4 . Fig- 
ure 1 shows how the 3PCF in the non-linear regime can be efficiently reconstructed from 
its lowest harmonic multipoles. In particular, for the isosceles triangle with side length 
r = 1 Mpc/h, only three harmonics contribute significantly to the 3PCF. The monopole 
essentially determines the amplitude, the (negative) quadrupole and, to a lesser extent, the 
octopole shape the configuration dependence (see top left and middle panels). Thus the full 
configuration dependence of the 3PCF can be reconstructed from these three harmonics (see 



4 the four integrals correspond to a mass function integral, Legendre transform of the bispectrum, and a 
double Bessel integral to get from Fourier to real space 
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top right panel) to an excellent approximation. Only collinear configurations, 6 = 0, n have 
(increasingly smaller) contributions from higher order terms. Similarly, the triangles with 
7*1 = 2 Mpc/h, and r 2 = 6 Mpc/h are shaped by the quadrupole and therefore, the config- 
uration dependence is encoded in this single harmonic to a large extent (see lower panels). 
These examples are typical for all non-degenerate triangles we have checked. 

As we shall show below (see §4), in general, only the lowest multipoles are non-zero and 
thus full configuration dependence of the non-linear cosmological 3PCF can be compressed 
in a few numbers. Thus using the multipole expansion greatly reduces the problem of 
constraining non-gaussianity from gravitational clustering. 

4. Comparison to Numerical Simulations 

We have used two sets of ACDM simulations from the public Virgo simulation archive 5 
with cosmological parameters f2 = 0.3, Qa = 0.7, h = 0.7, T = 0.21 and cr 8 = 0.9 and no 
baryons. The original Virgo simulation has 256 3 particles in a box-size of L = 239.5/i~ 1 Mpc, 
mass resolution of 6.86 • 10 10 /i _1 M Q and softening length L so f t = 25/i _1 kpc, and a larger 
box (VLS) simulation L = 479/i -1 Mpc that contains 512 3 particles, same mass resolution 
than the original Virgo simulation and L so f t = 30/i _1 kpc . These simulations have been 
gravitationally evolved using a P3M code (Macfarland et al. 1998; Couchman et al. 1995). 

Bispectrum is computed through FFT's on 256 3 grid-points with the method of Scocci- 
marro et al. (1998) (see their Appendix A). The 3PCF is measured with the fast algorithm 
based on multi-resolution KD-trees of Moore et al. (2001) and Gray et al. (2004) using the 
estimator of Szapudi & Szalay (1998). Limited by computational resources, we dilute simu- 
lations to 10 percent for the 3PCF at scales r ^ 2 Mpc/h and at 1 percent at larger scales. 
For these computations the number of points in auxiliary random sets are roughly 10 times 
larger than those of diluted simulation data sets. The large box simulation is cut into eight 
independent sub- volumes (octants) of half box size and measured separately. Mean values are 
obtained by averaging these eight sub-volumes and the orginal Virgo simulation. Errorbars 
are simply computed from the dispersion over sub- volumes. We note that this error estimate 
might underestimate true errors 6 since sub-volumes are somewhat correlated. However, for 
the small scales analyzed in this paper, our sub-samples are effectively independent. 



5 http://www. mpa-garching.mpg.de/Virgo 

6 Using a large number of very large box simulacions should yield a more correct error estimate on large 
scales 
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Our approach is to discuss halo model predictions for a fiducial model that provides a 
good fit to the 2-point statistics (i.e. power spectrum and 2PCF) measured in the simulations. 
Then we explore systematically the ability of such a model of describing 3-point statistics 
in N-body experiments. We set cosmological parameters to match those of our simulations 
and the halo concentration parameter is set to Co = 9 and (3 = 0.13 in eq(3), as in SSHJ. We 
have verified that changing the amplitude and slope of the concentration parameter within 
a reasonable range (i.e. by a factor ~ 2) has a negligible effect on the 2PCF and 3PCF 
for k <, 10, or r >, 0.2Mpc/h (see also Fig.2 in Takada & Jain 2003a), therefore we fix this 
parameter to its fiducial values to simplify further studies. 

Basic halo model parameters yield a 2PCF that does not fit well our N-body results. 
This is clearly seen in Figs. 2-6, especially for scales, r ss 1 Mpc/h. To improve the precision, 
we introduce standard "tweaks" of the basic set of ingredients of the model. The assumption 
that the halo boundary (the radius up to which halo mass is encompassed) is exactly given 
by the virial radius is in fact quite arbitrary. Treating this as an additional free parameter 
to the 1-halo terms in the 2PCF and 3PCF improves the transition scales, 0.5 — 5 Mpc/h, 
between linear (or quasi-linear) to fully non-linear scales (see Figs. 2 and 5). Fig. 5 shows that 
a halo boundary beyond the virial radius increases the 2PCF on small scales and reduces 
the "bump" in the reduced bispectrum for equilateral triangles Q eq (k). Setting the bound- 
ary to the fiducial value, 1.3R V: yields good agreement both with the 2PCF and 3PCF in 
simulations. In particular, for the 2PCF, the agreement with N-body is comparable to the 
fitting function provided by Smith et al.(2003). However, on large scales, halo model tends 
to slightly overpredict simulations, whereas Smith et al. slightly underpredicts them. Simi- 
larly, Fig. 4 shows that our fiducial model provides a good fit to the matter power spectrum 
measured in our N-body simulations. 

Our results are in full agreement with that of Takada & Jain (2003a) who also imple- 
mented the "exclusion effect" . This excludes correlations between haloes at distances smaller 
than the sum of their virial radii. Halo exclusion is an inherently real space phenomenon, 
and it would be difficult, if not impossible, to correctly implement in our Fourier-based ap- 
proach. We have checked using approximations (which essentially cut off the two-halo term 
on small scales) that halo exclusion plays a subdominant role in describing the 3PCF (see 
also Fig.8 in Takada & Jain 2003a). 

In addition, we try to mimic the finite volume effects affecting N-body simulations 
by imposing a cut-off in the mass function, excluding haloes of masses M £ M cut . This 
has a critical effect on clustering on scales r £ 5Mpc/h. We find that imposing no cut-off 
overpredicts the clustering observed in simulations whereas excluding haloes of M £ 1O 14 5 M 
results in too low a 2PCF. This is even more evident in the reduced bispectrum for equilateral 
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triangles Q eq (k) (see Fig. 5). Our fiducial value M cut = 1O 15 M reproduces reasonably well 
simulation results. This behavior was already observed by SSHJ and Wang et al. (2004) 
who found different best-fit values for their smaller-box (L = 300Mpc/h) simulations. While 
these works did not discuss the effect of the mass cut-off on the 2PCF, we have made sure 
that our fiducial choice, motivated by the 3PCF, still reproduces the 2PCF accurately. 

The above two extra parameters fix our fiducial model. In what follows, we shall use it 
to work out detailed predictions for the 3PCF and its harmonic multipoles Q (see eqs(ll) 
and (12)), and compare them to measurements N-body simulations. Fig. 7 displays the 
reduced bispectrum for several different triangular configurations. It appears that large 
scales present a more pronounced configuration dependence than predicted by the three- 
halo term. Triangles on small-scales are rather flat, except for the "shoulder" displayed by 
the isosceles triangles. 

Halo model prediction for the reduced 3PCF for equilateral triangles, Q eq (j') is in good 
agreement with N-body results as shown in Fig. 8. On scales r £ 3Mpc/h the model slightly 
overpredicts simulations, similarly to the 2PCF on the same scales (see Fig. 2 and Fig. 3). 
Note that the non-linear halo model converges, albeit slowly, to the PT on large-scales. 

Fig.9 displays the configuration dependence of the 3PCF for the transition scales as 
predicted by the halo model (lines) compared to numerical simulations (symbols). The 
agreement between the model and simulations is generally within the errorbars. Moreover, 
the rich configuration and scale dependence made up by the non-trivial interplay between 
different halo terms is observed in simulations. This can be considered as the best confir- 
mation that the basic idea of halo models carries over to the three-point function without 
major modifications. 

Fig.10 shows the reduced 3PCF Q(r) for the same triangles whereas Fig. 11 shows an 
interpolated surface displaying the (9, q) dependence for the mean values of the N-body for 
r = 2 Mpc/h. Note that since Q(r) is predicted to be of order unity for all configurations, 
it is more suitable for high-precision comparisons on linear scale plots. On the down side, 
the complex structure of this ratio statistic combining 2PCF and 3PCF complicates inter- 
pretation. In particular, halo model describes the configuration dependence of the reduced 
3PCF accurately for r £ 0.5Mpc/h, but it becomes increasingly inaccurate for larger scales. 
The slight overall amplitude mismatch observed (see e.g., the case for isosceles triangles, 
q — 1, in left column of Fig.10) is due to the halo model 2PCF overpredicting simulations. 
We have checked that using fitting formula of (Smith et al. 2003) for the 2PCF instead of 
the halo model does not improve the agreement with simulations significantly. Note that for 
intermediate angles theoretical predictions always agree with the simulations, in line with 
our earlier results for equilateral triangles (see Fig. 8). The most visible difference between 
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model and N-body results is the lack of configuration dependence predicted by the halo 
model on larger scales, r >, 0.5Mpc/h. 

The multipole expansion of the 3PCF, eq.(12), is particularly useful to pinpoint why 
the halo model fails to reproduce numerical results. Figs. 12-15 show the 3PCF multipoles 
for different triangles and scales r = 0.25 — 2 Mpc/h. In general, the amplitude of the 
coefficients falls off rapidly with multipole order. This makes multipoles especially convenient 
to compress information in the cosmological 3PCF in the non-linear regime (see also Szapudi 
2004). 

Isosceles triangles are exceptional: on small angles the third side of the triangle is very 
small, which results in a rapid increase. This can be described in multipole space only 
with a high number of multipoles (such as a Dirac 5-function would have infinite number of 
multipoles). The situation worsens towards larger scales, due to even more rapid increase 
of the 3PCF on small angles. While this technical deficiency can be overcome by simply 
smoothing the correlation function over angles ("band limiting"), we note that on scales 
above r >, 5Mpc/h, where this effect becomes severe, PT becomes more and more accurate, 
therefore halo models are not necessary for dark matter predictions (Fry 1984; Jing & Boerner 
1997; Barriga & Gaztahaga 2002). 

More interestingly, for r £ 0.5Mpc/h, halo model consistently underpredicts the quadrupole 
moment (and to a lesser extent, the octopole). The observed lack of configuration dependence 
in the halo model with respect to simulations is primarily due to a quadrupole deficit in the 
prediction (see Figs.9 and 10). This might be caused by the 2-halo term, which shows hardly 
any variation with angle on large-scales r £ 0.5Mpc/h (see Figs.9 and 10), contrary to what 
one would expect for halo-halo spatial correlations. One possible reason for this quadrupole 
deficit in the two-halo term is that our implementation of the halo model assumes that haloes 
are spherical, whereas real virialized objects in high-resolution CDM simulations do have a 
typically asymmetric (triaxial) shape (see e.g., Barnes & Efstathiou 1987; Jing & Suto 2002; 
Moore et al. 2004). SSHJ already pointed out in their bispectrum analysis that relaxing the 
sphericity hypothesis for halo shapes could bring halo model to a closer agreement with the 
configuration dependence observed in numerical simulations. Although our analysis suggests 
a similar conclusion, it is unclear whether one-halo (as pointed out by SSHJ) or rather two- 
halo terms (that show a significant lack of configuration dependence at all scales) should 
carry the missing quadrupole. 

On the other hand halo substructure (Sheth 2003; Sheth & Jain 2003) could play a 
significant role on scales comparable to large-cluster sized haloes IMpc/h. However, a recent 
analysis (Dolney et al. 2004) suggests that substructure tends to attenuate the amplitude of 
the reduced bispectrum on small scales. This would render our analytic predictions, at least 
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our fiducial model, in disagreement with N-body measurements (see Figs. 5 and 6). 

Clearly, there are other potential improvements to our implementation of the halo model, 
most importantly a consistent definition of the mass function for a modified halo boundary 
(see White 2002; Hu & Kravtsov 2003), or allowing for a steeper inner halo profile (see Moore 
et al. 1998) what could bring model predictions closer to measurements in high-resolution 
N-body simulations. We plan to address these issues in future work. 

5. Conclusions 

The revival of the halo model in recent years has been triggered by the ability of high- 
resolution N-body simulations to test theoretical predictions with great precision. In this 
paper we have presented a detailed comparison of halo model predictions against simulations 
for the 2PCF and 3PCF in real and Fourier space. Our analysis has focused on the transition 
scales, r « 0.25 — 5Mpc/h that connect large (quasi-linear) scales, appropriately described by 
PT, with highly non-linear scales where phenomena such as the stable clustering hypothesis 
require higher resolution simulations to be accurately tested. 

Our results show that halo boundary and mass function cut-off have a significant effect 
on the three-point correlation functions on these transition scales, and thus these statistics 
can be used to constrain such halo model parameters. A fiducial model with halo boundary 
r = 1.3R V Mpc/h and mass cut-off M cut = 10 15 M Q brings theoretical predictions in close 
agreement to our simulations. The success of our fiducial model in explaining non-linear 
gravitational clustering has been comprehensively demonstrated in Fourier space for the 
power spectrum (Fig. 4), the reduced bispectrum (Figs. 5 and 6), as well as in real space, 
through the two-point correlation function (Figs. 2 and 3), 3PCF (Fig. 9), and the reduced 
3PCF (Figs. 8, and 10) along with its multipole moments (Figs. 12-15). 

Halo model predictions are in good agreement with the configuration and scale de- 
pendence of the reduced bispectrum at all the scales measured in N-body results, k = 
0.25 — 1 h/Mpc. On the other hand, while the model correctly predicts the amplitudes of 
the 3PCF on intermediate angles for all the scales tested, it exhibits a lack of configuration 
dependence that becomes more significant on scales larger than the largest halos seen in 
simulations r £ 1 Mpc/h. Although the reason for this is not clear, we suggest that non- 
spherical haloes could produce a more pronounced configuration dependence in theoretical 
predictions. This could be realized through the two-halo terms that show a rather flat be- 
havior in the current implementation of the model. However if the orientation of the two 
halos is random, the non-sphericity should cancel out. This issue certainly deserves further 
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attention in future analyses of the dark-matter and galaxy N-point correlation functions. 

It is particularly useful to decompose the reduced 3PCF in its harmonic multipoles Qes, 
as given by eq(17). We have seen that only the lowest orders have non-zero amplitude and 
thus the non-linear 3PCF can be compressed in a few numbers (see Fig. 1). Furthermore, 
assuming a model for the halo occupation distribution and an appropriate modeling of red- 
shift distortions, one should be able to formally decompose the 3PCF of galaxies in the 
exact same way. In particular, observables such as the harmonic multipoles of the galaxy 
3PCF measured in large volume surveys should be largely uncorrelated and approximately 
Gaussian distributed, in a similar way it happens for the angular power spectrum multipoles 
of the CMB anisotropy Cg. 

Thus the Q/s must be tightly related to the signal-to-noise eigenmodes of the 3PCF 
(Scoccimarro 2000; Gaztanaga & Scoccimarro 2005). Moreover Gaztanaga & Scoccimarro 
(2005) have shown that the highest signal-to-noise Q-eigenmodes measured in mock galaxy 
catalogs can be efficiently used to recover bias parameters. It is remarkable that the first three 
Q-eigenmodes have a configuration dependence that is extremely similar to the monopole, 
dipole and quadrupole terms (compare their Fig. 10 with middle panels in Fig.l of this paper). 
This suggests that the Q^s are indeed closely related to the uncorrelated modes of the galaxy 
3PCF and thus they may simplify the procedure of constraining galaxy bias from the 3PCF. 

Alternatively, the proposed multipole approach can be readily applied to other clustering 
statistics of the mass density field in the non-linear regime such as the projected mass 3PCF, 
that is barely affected by redshift distortions (see Zheng 2004, for a recent implementation 
using Fourier series), or the 3PCF of the convergence field that probes the lensing potential 
(see e.g., Takada & Jain 2003b). 
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Fig. 1. — Reconstructing the 3PCF form its harmonics. Only the lowest harmonic coefficients 
are non-zero (see left panels). Middle panels show the configuration dependence of individual 
harmonic terms i.e. harmonic coefficients times the Legendre polynomials. We display terms 
I = 0, 1, 2 and 3 using increasingly thicker lines. Right panels show the sum over harmonics 
up to a given £ ma x = 1,2,3 and 10 (from thin to thick lines). In general, the non-linear 
3PCF is dominated by the monopole, that sets the amplitude, and quadrupole, that shapes 
the configuration dependence. Top panels show isosceles triangles of length r = 1 Mpc/h, 
while lower panels display more elongated triangles with basic length r = 2 Mpc/h and 
length-ratio q = 3. 
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r [Mpc/h] 

Fig. 2. — Halo boundary effects on the 2PCF: defining halo boundary at r = R v underes- 
timates £(r) on the transition scales. Larger boundaries yield a larger one-halo (lh) con- 
tribution. For r = 1.3 R v (solid line), we find good agreement with N-body results. The 
inset shows a blow-up of the transitions scales, where our fiducial model agrees well with 
numerical results and the Smith et al. (2003) fitting function. On larger scales the fiducial 
model overpredicts numerical simulations, but converges to them for r ^ 20Mpc/h. 
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Fig. 3. — Same as Fig.2 for the mass cut-off: in order to mimic the absence of haloes of mass 
M > M cut in simulations, our fiducial model incorporates a cut-off in the mass function, 
M cut = 1O 15 M (solid line), what lowers the 2PCF on non-linear scales (i.e. those dominated 
by the one-halo term). The Smith et al. fitting function is only shown in the inset for clarity. 
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k [h/Mpc] 

Fig. 4. — Mass power-spectrum per logarithmic wavenumber, A 2 (k) = k 3 P(k)/2n 2 , in the 
fiducial halo model (solid line, see text for details) compared to N-body simulations (sym- 
bols). One (lh) and two (2h) halo contributions are shown to have comparable amplitude 
in the transition scales k « 0.5 h/Mpc. For reference, it is also shown a fitting function to 
numerical simulations (dashed line; Smith et al. (2003)). Halo model agrees very well with 
N-body except on scales k « 0.1 — 0.5 h/Mpc, where it tends to slightly overpredict numerical 
simulations. Note that we use a box-size L = 239.5 Mpc/h for small scales, k ^ 2 h/Mpc, 
whereas a larger box (VLS) simulation L = 479 Mpc/h is used for larger scales, k £j 2 h/Mpc. 
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Fig. 5. — Halo boundary effect on the reduced bispectrum for equilateral triangles Q eq {r). 
Qeqir) is very sensitive to changes in the halo boundary: larger boundaries enhance more the 
denominator (that depends on the 2PCF) than the 3PCF, what results in a significant drop 
of the reduced bispectrum on non-linear scales. On large-scales PT holds (see dotted line). 
Note that the fiducial model matches simulations, but exhibits a bump on k ~ 1 — 2h/Mpc, 
that is not seen in N-body results. 
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Fig. 6. — Same as Fig.5, but for the mass cut-off. Excluding the largest mass haloes signifi- 
cantly lowers the amplitude of the reduced bispectrum. 
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Fig. 7. — Reduced bispectrum for different triangle configurations. Rows show increasingly 
larger scales from top to bottom, while columns evolve towards more symmetric triangles 
from left (q = 4) to right (q = 1, isosceles triangle), where q = k 2 /ki, and k = k\. Lines 
show the total (solid), one-halo (short-dashed), two-halo (long-dashed), and three-halo (dot- 
dashed) terms. Note the continuity of all the lines through plotted triangle configurations 
(from left to right). This "boundary condition" can be expressed as Q(k,q,6/-K = 1) = 
Q(k,q — 1,9/tc = 0), since these configurations describe the same triangle: the fact that 
the lines continue as expected shows the accuracy of our integration. It is seen that large 
scales show a more pronounced configuration dependence as dictated by the three-halo term. 
Triangles on small-scales are rather flat, except for the "shoulder" displayed by the isosceles 
triangles. 
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Fig. 8. — Reduced 3PCF for equilateral triangles: the fiducial halo model (upper solid line) 
shows good agreement with simulations. The small-scale enhancement is determined by the 
one-halo term, while the amplitude on the transition scales r « 1 — 3 Mpc/h reflects the non- 
trivial interplay between different halo terms (lower solid lines). The good agreement found 
with N-body in this range is a remarkable success of the halo model. The observed mismatch 
between model and simulations on larger scales is due to the model slightly overestimating 
the 2PCF (that goes into the denominator of Q eq ) on that range of scales (see Fig.2). 
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Fig. 9. — 3PCF for different triangle configurations: lines go as in Fig. 7. Halo model predic- 
tions are shown to be in good agreement with simulations (symbols). On scales r £ 1 Mpc/h 
the model shows a flatter configuration dependence than simulations, but they still agree 
within errorbars for most of the cases. Note that in real space the triangle "boundary con- 
dition" (see text and Fig.7 for details) is expressed as ((r, q, 9/tt = 1) = £(r, q + 1, 9/n = 0), 
where q = r 2 /r 1; and r = r\. 
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Fig. 10. — Same as Fig. 9 but for the reduced 3PCF. Halo model reproduces simulations 
reasonably well, specially for intermediate angles, what is in agreement with Fig.8 that 
shows equilateral triangles alone. The discrepancy with simulations becomes more evident 
in this ratio statistic, as it combines the complicated configuration dependences of both 
2PCF and 3PCF. As in Fig.9, the discrepancy becomes more significant for r ^ 0.5 Mpc/h. 
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Fig. 11. — Smoothly interpolated surface showing the configuration dependence of the re- 
duced 3PCF at r = 2 Mpc/h from N-body simulations. It exhibits a symmetric convex 
shape, what is in qualitative agreement with predictions from the three-halo term in the 
halo model. For isosceles triangles this symmetry is broken at small-angles due to non-linear 
effects, as reproduced by the one-halo contribution (see Fig. 10). 
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Fig. 12. — Harmonic multipoles of the 3PCF at r = 0.25 Mpc/h: theoretical predictions 
are in excellent agreement with simulations. Only lower order multipoles (£ £ 5, except for 
isosceles triangles) have a non-zero contribution. 
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Fig. 13. — Same as Fig. 12 but for r = 0.5 Mpc/h. 
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Fig. 14. — Same as Fig. 12 but for r = 1 Mpc/h. Note how the higher multipoles becomes sig- 
nificant for isosceles triangles (q = 1), and that theoretical predictions tend to underestimate 
the quadrupole moment for q > 1. 
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Fig. 15. — Same as Fig. 12 but for r = 2 Mpc/h. In line with Fig. 14, it shows that halo 
model clearly underpredicts the quadrupole (and to a lesser extent, the octopole) observed 
in numerical simulations. This quadrupole deficit explains the lack of ^-dependence in halo 
model predictions for Q(r, q, 8), as shown in the lower rows of Fig. 10. 



